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SUMMARY 


The spatial distribution of the numerical disturbances that are generated dur- 
ing a numerical solution of a flow, is examined. It is shown that the distribution of 
the disturbances is not uniform. In regions where the structure of a flow is simple, 
the magnitude of the generated disturbances is small and their decay is fast. How- 
ever in complex-flow regions, as in separation and vortical areas, large- magnitude 
disturbances appear and their decay may be very slow. The observed nonuniformity 
of the numerical disturbances makes possible the reduction of the calculation time 
by application of, what may be called, the partial-grid calculation technique, in 
which a major part of the calculation procedure is applied in selective subregions, 
where the velocity disturbances are large, and not within the whole grid. This 
technique is expected to prove beneficial in large-scale calculations such as the flow 
about complete aircraft configurations at high angle of attack. Also, it has been 
shown that if the Navier-Stokes equations are written in a generalized coordinate 
system, then in regions in which the grid is fine, such as near solid boundaries, the 
norms become infinitesimally small, because in these regions the Jacobian has very 
large values. Thus, the norms, unless they are unsealed by the Jacobian, reflect 
only the changes that happen at the outer boundaries of the computation domain, 
where the value of the Jacobian approaches unity, and not in the whole flow field. 


1. INTRODUCTION 

For the numerical calculation of a flow, a reasonable initial solution is selected 
and an iterative process is applied until, by means of one or more convergence 
criteria, a converged solution is reached. The iterative procedure may last a few 
steps, or hundreds, or even thousands of steps, depending mainly on the nature 
of the flow, but also on the fineness of the grid and on the type of the numerical 
scheme. 

The convergence criteria that are currently used are numerical or physical ones. 
Thus, during a numerical solution, one may observe the variation of one of the norms 
of the residual of the finite- difference equations, or the variation of the difference 
of the calculated quantities, from step to step, or, if appropriate, one may check 
the value of a physical quantity, such as the lift or the pressure coefficient. In the 
first case, the calculation terminates when the norms or the differences reach an 
assumed “zero” value, while, in the latter case, when no more change is observed 
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in the value of the coefficient. 


The classical convergence criteria reflect the variation of the mean value of a 
quantity and they do not give information about the distribution of this quantity 
within the examined flow field. However, the flows in general are nonuniform. 
Especially large nonuniformities are observed in separated and vortex flows, in flows 
near boundaries, in shear flows, and in shock waves. In all these cases the kinematic 
and, probably, the thermodynamic quantities are subjected to a large variation, 
often in a region that is just a small portion of the whole flow field. Therefore it 
is expected that during the development of a numerical solution more steps will be 
required for the decay of the disturbances in these “complex” flow-regions, rather 
than in the ones which have a simpler structure. 

In this paper it will be shown that indeed, the intensity of the generated dis- 
turbances to the calculated flow quantities, and the time required for their decay 
is much higher in the complex-flow regions than in the simple-structure regions. 
Also, it will be shown that instead of the standard repetition of the numerical inte- 
gration across all the extent of the computational domain (grid), until the solution 
converges, it is possible to do a part of the calculation cycle selectively, in subre- 
gions characterized by large values of the disturbance field. In this case the overall 
time required for the calculation of a flow will be smaller. Considering the present 
requirements for large-scale calculations about bodies at high incidence, this possi- 
bility is attractive. That because these types of flows are characterized by a rather 
simple structure in the windward side, and by a complex structure in the leeward- 
side. Also, the computer time required for the estimation of the flow field about a 
modern air vehicle is of the order of tens of hours in a supercomputer. Thus there 
is a need to reduce the required computer time. 

For the calculation of the flows which are used for the study of the numerical 
disturbance fields, the code F3D of NASA Ames (ref. 1) was used. The code F3D 
is a two-factor, implicit, approximately factored, finite- difference code which solves 
a conservative form of the thin-layer Navier-Stokes equations cast in generalized 
coordinates, so that it can be readily used for computing flows about complex 
configurations. 

In the generalized form, all the terms of the Navier-Stokes equations are divided 
by the transformation Jacobian of the system of coordinates. The value of the 
Jacobian within a grid may vary by orders of magnitude. Especially large values 
are observed near walls, where usually the grid is fine, while at the outer layers 
of a grid the Jacobian tends to unity. It will be shown that this variation may 
have a significant masking effect on the numerical mean- value convergence criteria, 
if the various quantities on which these criteria are based are not unsealed by the 
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Jacobian. By combining the mean- value convergence criteria with surveys of the 
disturbance fields it was found that, unless the criteria are unsealed, the level of the 
disturbances near the wall is not shown. 


2. GOVERNING EQUATIONS AND CALCULATION PROCEDURES 

The conservation equations of mass, momentum and energy can be represented 
in a flux- vector form that is convenient for numerical simulation as (ref. 1): 

drQ + d((F + F v ) + d n (G + G v ) + d c {H + F«) = 0 (1) 

where r is the time, and the independent spatial variables £,77, and £ are chosen 
to map a curvilinear, body-conforming discretization into a uniform computational 
space. In the system used in this study, ( denotes the curvilinear axis in the direction 
of the main body axis, 77 denotes the circumferential and ( thejiormal to thejbody 
direction. As opposed to the inviscid flux terms F, G, and if, the terms F V ,G V 
and H v are fluxes containing the viscous derivatives. A nondimensional form of the 
equations is used throughout this work. Lengths are scaled by the length of the 
body, L, velocity components by the free-stream velocity of sound, aoo> the pressure 
and the total energy per unit volume (e) by />oo«oo 2 * The other quantities, such as 
T, pi fi are scaled by their free-stream values. 

For body-conforming coordinates and high Reynolds number flow the thin-layer 
approximation can be applied (ref. 2, 3, 4). 

d T Q + dtF + diG + dcH = Re-'d^S (2) 

where only viscous terms in the normal to the body direction are retained. These 
have been collected into the vector S and the nondimensional Reynolds number Re 
is factored from the viscous flux term. 

An implicit approximately factored scheme for the thin-layer, Navier-Stokes 
equations that uses central differencing in the 77 and £ directions and upwind differ- 
encing in the £ direction can be written in the form: 

[/ -I- hS b ( (A + ) n + h8<C n - hRe~ 1 8^J~ 1 M n J - Di\ ( ] 
x[I + hS((A~) n + h8r,B n - £>i|„]AQ n = 

- At{S b [(F+r - F+] + */[(F-)" - F-] (3) 

+ 6 v (G n - Goo) + k(H n - £00) - - S^)} 

- D e (Q n - § 00 ) 
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where h = At or At/2 for first or second order time accuracy, and the free-stream 
base solution is used. Here 6 is typically a three-point, second-order-accurate, cen- 
tral difference operator, 6 is a midpoint operator used with the viscous terms, and 
the operators 8 | and 8^ are backward and forward three-point difference operators. 

A. A A 

The flux F has been split into F + and F~ , according to its eigenvalues(ref. 5), and 
the matrices A ± ,B,C, and M result from local linearization of the fluxes about 
the previous time level. In equation (3), J denotes the Jacobian of the coordinate 
transformation. Dissipation operators, D e and D{, are used in the central space 
differencing directions(ref. 1). 


The factored left-hand-side operators can be readily solved by sweeping in the £ 
direction and inverting tridiagonal matrices with five by five blocks. This two-factor 
implicit scheme is readily vectorized in planes of £ = constant. 


Concerning the calculation procedure, during a numerical solution of a flow, it 
is known that in case of nonsteady flows a time-accurate method must be employed. 
After the selection of some reasonable initial solution, an integration in time is done, 
with time steps commensurate with the nonsteady phenomenon that is calculated. 
In steady state calculations, usually non-time-accurate techniques are applied, that 
bring the flow from an arbitrary initial guess to its asymptotic solution, with the 
least possible computational work. The use of a space varying time step is a classical 
technique, in this sense. 


Often, for simplicity, the impulsive start, i.e. the sudden immersion of a body in 
a uniform flow, is used. Various techniques have been developed that accelerate the 
convergence. For both steady and nonsteady flows, a faster convergence is reached 
if instead of the impulsive start, a solution of a similar flow is applied initially. Also, 
in the case of a nonsteady flow, it is better to calculate first an approximate steady 
solution and then to apply the time-accurate calculation. One way to accelerate 
convergence to a steady state is to obtain a good initial guess for a fine mesh by 
first iterating on a sequence of coarse grids and then interpolate the solution up to 
the next refined grid (see, for example, ref. 6). This is called the mesh sequence 
procedure. 


Some of the described techniques of accelaration of the convergence will be 
applied in the numerical solution of the flows that are used in this work, for studying 
their effect in the numerical disturbance field. 
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3. DESCRIPTION OF THE CALCULATED FLOWS 


In this section, the flows that have been used for studying the spatial varia- 
tion of the disturbance field will be briefly described. The detailed results of the 
investigated flows will be reported in other particular publications. 

The first case is related to the calculation of the turbulent flow about a prolate 
spheroid, at a = 30°, Moo = 0.25 and Rei = 43x10®. The prolate spheroid is a 
body of simple geometry that can provide significant understanding of the complex 
flow that can develop about axisymmetric bodies at incidence. In a cross-section 
of such a flow, under the action of circumferential pressure gradient, the outer flow 
approaching the windward plane of symmetry turns and flows outwards along the 
body, from the windward toward the leeward side. The boundary layer which is 
formed in this way, separates from the body at a point on the leeward side. The fluid 
then leaves the body along a feeding sheet and rolls up to form a primary vortex 
system on the leeward side of the body. The pair of the primary vortices induces 
a flow toward the body surface, which at the point of attachment turns outwards 
toward the windward side. For primary vortices of sufficient strength, a secondary 
separation is induced to the boundary layer below them. Thus, a secondary vortex 
structure may appear that rotates in an opposite direction than the primary one. 
Both the vortex systems are carried downstream by the axial component of the 
flow. 


Results of the calculations are shown in figure 1. The primary vortex structure 
is clearly shown in the velocity-vectors plot (fig. lb), while the existence of the 
secondary- vortex is detected more easily in figure la, where surface particle traces 
are used to simulate the surface shear-stress lines visualized in oil-flow experiments. 
These lines converge toward the primary and the secondary separation lines and 
diverge from the attachment line, which lies between them. 

The calculation was non-time-accurate, with the time step varying in space, 
as a function of the local Jacobian. No similar solution was available, so the flow 
was started by imposing impulsive conditions. The grid that was used consisted of 
121 points along the direction of the main axis of the body, 100 points around its 
circular section and 65 points radially from the surface. Initially, for accelerating the 
solution, a coarse grid was used. That grid was formed by considering each second 
point of the regular grid, in the streamwise and in the circumferential direction. 
Also, the values of the coefficients of the numerical viscosity terms of equation (3), 
initially were rather large, also for the purpose of acceleration of the convergence. 

A final change that was imposed during the numerical solution was related to an 
improvement of the turbulence model. The standard model used in the code F3D, is 
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the two-layer, Cebeci type, algebraic eddy-viscosity model reported by Baldwin and 
Lomax (ref .2). In this model the value of the turbulence viscosity coefficient, fi t , 
at a particular velocity profile is a function of the maximum of the product of the 
vorticity and the distance from the wall. Degani and Schiff (ref. 7) have observed 
that in flows which include crossflow vortices there exist two maxima in the distri- 
bution of the vorticity product in the velocity profiles, one caused by the boundary 
layer and one caused by the overlying vortical structure. The consideration of only 
the lower maximum gives better results. In the calculation of the flow about the 
prolate spheroid, initially no distinction had been made between the maximum of 
the vorticity product caused by the boundary layer and the one caused by the vor- 
tex. Later in the calculation, this distinction was considered. This change affected 
strongly the separated vortical region. The primary vortices became larger and the 
secondary vortices, which previously had not shown, appeared. The attached- flow 
was not affected by this change because there the value of the fit was kept the same. 


All of the forementioned changes that were applied during the calculation had 
an impact in the state of the numerical solution, in the sense that when one of these 
changes was imposed, the solution was perturbed. Thus, an excellent opportunity 
appeared for studying the spatial growth and decay of the velocity perturbations 
during the process of convergence. 


The second case that is used in this work, is related to the impingement of a 
jet, in the presence of a crossflow. In the simulation, the jet was originated from a 
nozzle and it was directed perpendicularly to the ground. The ratio of the crossflow 
velocity to the mean jet velocity was equal to 0.223, and the Reynolds number, 
based on the diameter of the nozzle, was 4.4x10 s . The flow was calculated as a 
turbulent one. Van Dalsem (ref. 8) has investigated the case where the nozzle has a 
circular cross-section. The present author has studied the case of a jet of elliptical 
cross-section (ref. 9). The axes ratio has been variable and the major axis has been 
aligned or normal to the crossflow. In this paper, the case of axes ratio equal to 
1.5, with the major axis aligned with the crossflow, will be used for the study of 
the numerical disturbance field. 


In this case the calculations were time accurate and the solution of the round 
jet was used as initial one. After changing the grid, to account for the new cross- 
section, the calculation started with the assumption that at each point the flow 
quantities were those of the round jet. The resulting perturbation and its gradual 
decay will be analysed in the next section. To give an idea of the particular flow, in 
figure 2a, the trajectories of particle traces released from the nozzle face are shown. 
Also, in figure 2b the velocity vectors are shown in the vertical symmetry plane. 
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4. RESULTS 


4.1 Generation and Decay of the Disturbances 

For the study of the generation and decay of the disturbance field during the 
computation of the flows that have been selected, mean- value, as well as distribution 
convergence criteria will be used. The It norm, or Euclidean norm, is the most 
widely used mean- value criterion. It reflects the change of the residual of the finite- 
difference equations. In this paper, the I 2 norm will be used, scaled or unsealed 
with the Jacobian. In addition, the mean- value of the A Q vector will be used. 

The change of the velocity field (Atii,i = 1,2,3) between two successive steps, 
has been selected for studying the spatial variation of the disturbance of a numerical 
solution. The A m is defined by the relation u< n+1 = Ui n + Au{. In what follows, 
A Ui will be called the velocity disturbance. When the velocity disturbance during 
a numerical solution of a flow is equal to zero, it is assumed that the flow solution 
has converged. In the various velocity plots that will be presented in this paper, 
the At n have been magnified by a factor of 10 3 -10 5 , so that their change will be 
visible. In this type of plot the vectors may cross a solid boundary. 

Starting with the case of the prolate spheroid, the history of the generation 
and decay of the mean value of the AQ-disturbances is shown in figure 3. It is 
seen in figure 3 that the initial strong perturbation, because of the assumption of 
uniform flow conditions gradually decays, but not smoothly. Three large discontin- 
uous increases are observed in the velocity disturbance. These discontinuities reflect 
the perturbation of the numerical solution each time that a change was imposed. 
The first discontinuity appears at N=250 (number of iterations), where the coarse- 
grid solution was interpolated into the fine grid. The next change was imposed 
at N=530, and it was related to the reduction of the values of the coefficients of 
the numerical dissipation terms. Again, this disturbance has shown up in the A Q 
curve, by a considerable jump. However the largest disturbance appeared when the 
turbulence model changed at N=1830. It is seen in figure 3 that the A Q curve 
presents a discontinuous increase of two orders of magnitude at this point in the 
convergence history. 

Concerning the development in time of the spatial variation of the numerical 
disturbance field, the present results suggest that when a perturbation is imposed 
in the solution, large disturbances are generated at those points of the flow (grid) on 
which the new values imposed are different from the previous ones. During the im- 
mediate subsequent cycles of calculation, the disturbances are propagated outwards 


7 


from the initial regions of disturbance and gradually cover a large volume, which is 
not necessarily equal to the whole calculation domain. Then the disturbances start 
to decay, so that gradually the perturbation domain shrinks towards the initial per- 
turbation region and from there towards low value, in case of a steady flow. The 
present data has not covered the case of the nonsteady flows, however it seems that 
in that case numerical disturbances will always be present in the unsteady part of 
the flow, but not in the whole flow field. 

Of great importance for practical applications is the fact that the regions in 
which the numerical disturbances are large coincide with the regions in which the 
flow has large gradients. Thus it is easy to anticipate where it will take more time 
for a numerical solution to converge. 

A clear demonstration of the nonuniform generation and decay of the numerical 
disturbances is given by figure 4, in which the history of the development of the 
numerical disturbance field is given, in the case of the reduction of the value of the 
numerical dissipation terms. It is seen in figure 4 that the disturbances initially 
appear near the surface of the body (fig. 4a). Then they propagate and cover the 
boundary layer and the vortical separation region at the leeward region (fig. 4b); 
in the rest domain, “silence” conditions prevail. Gradually the disturbances die-out 
and uniform conditions of “silence” are established everywhere (figs. 4c-e). 

As an additional evidence of the nonuniform convergence of a numerical solu- 
tion, in figure 5 the impact of the modification of the turbulence model, at the time 
step N=1830, on the numerical disturbance field is shown. It has been mentioned in 
section 3 that this modification has affected only the separated-flow region and not 
the attached one. In figure 5 the same is seen to happen to the numerical distur- 
bance field. The disturbances appear only in the leeward vortical region, where the 
crossflow separation exists, propagate outwards to account for the enlargement of 
the vortices, and gradually decay when equilibrium conditions have been reached. 
In the windward region, where the flow is attached, the numerical solution is not 
disturbed at all. 

In the case of the impinging jet in the presence of a crossflow, the behaviour 
of the numerical disturbance field was found to be exactly similar with that of the 
flow about the prolate spheroid. The evolution of the disturbance field is shown in 
figure 6. As soon as the geometry of the grid is perturbed, by changing the circular 
nozzle into one of elliptic cross-section, large disturbances appear on the upstream 
face of the nozzle and, mainly, in the region of the jet (fig. 6a). The disturbances 
propagate radially but they do not cover the whole calculation domain (fig. 6b). 
The equilibrium conditions start to be established initially on the outer layers and 
gradually reach the jet (fig. 6c). After 600 steps (fig. 6d) the disturbance field 
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has been restricted to the vicinity of the jet and of the ground vortex and to the 
downstream wall flow. From then on, 1100 steps were required for the decay of these 
localized disturbances (fig. 6e). The final decay of the disturbances is an indication 
that the flow is rather steady. Also, Van Dalsem (ref. 8) has concluded that the 
round jet configuration, which he studied, shows only a small unsteadiness. Van 
Dalsem tried to develop a self-sustained unsteady jet in ground effect, by varying 
periodically the jet-exit velocity. However, once the forced pulsing was removed the 
flow quickly returned to a relatively steady-state condition. 


4.2. Effect of the Value of the Jacobian 

In the expression of the Navier-Stokes equations in a generalized system of 
coordinates all the terms are divided by the transformation Jacobian of the system 
of coordinates. However, the value of the Jacobian within a grid may vary by 
orders of magnitude. This variation affects the value of the numerical convergence 
criteria, because the disturbance quantities, on which the criteria are based, have 
been divided at the various points of the flow field by the variable Jacobian. Thus, 
the level of the disturbances in regions where the Jacobian has large values will be 
masked. This nondesired effect of the Jacobian will be examined in this section. 

The variation of the Jacobian in the crossflow and in the streamwise direction 
of the prolate spheroid is shown in figure 7. It it observed that the greatest variation 
of the Jacobian happens in the normal to the body direction. Thus, on the surface 
of the ellipsoid the Jacobian may be equal to 10 10 , depending on the axial position, 
while at the outer boundaries of the grid it may be even smaller than unity. In 
the streamwise direction the variation is somewhat smaller, the extreme values vary 
only three orders of magnitude. 

In the case of the jet in a cross-flow configuration high values of the Jacobian 
appear in the vicinity of the jet and near the impingement surface. The value of the 
Jacobian there is five to seven orders of magnitude greater than in the far field (fig. 
8). However, it is the jet and especially its impingement region where the greater 
values of the numerical disturbances occur (fig. 6). 

The masking effect of the Jacobian is very clearly shown if the values of the I 2 
norm are estimated, during a numerical calculation, with the residual of the Navier- 
Stokes equations scaled or unsealed by J. The results of this type of calculation are 
shown in figure 9, for the case of the prolate spheroid. The unsealed I 2 norm curve 
(fig. 9a) is similar in shape with the Au< curve (fig. 3). All the changes that were 
applied during the numerical solution are shown in this curve, by discontinuous 
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increases of the disturbance value. On the contrary, it is hard to detect the effect 
of the changes on the convergence level in the case of the standard Z 2 norm curve 
(fig. 9b). It is only the change of the turbulence model that somehow is detected, 
by a small discontinuous change at N=1830, and then by a gradual increase. It is 
reminded that in this particular application the disturbances start to grow near the 
surface and then propagate outwards (section 4.1). In the case of the change of the 
turbulence model, the outward extent was great because the vortex structure was 
strongly affected. 


The same observations apply, also, in the case of the jet in crossflow. The 
variation of the two I 2 norms is given in figure 10. It is seen that while the unsealed 
I 2 norm curve (fig. 10a) presents repeated, if not periodic, fluctuations, the standard 
I 2 norm (fig. 10b) is very smooth. An examination of the velocity-disturbance 
plots indicated that these fluctuations correspond to disturbance burstings within 
the jet at the impingement region near the wall. Most probably these burstings 
are of numerical nature and do not indicate unsteadiness of the flow. A detailed 
examination is necessary. It is evident from this example that the standard Z 2 norm 
is not appropriate to detect disturbances near the wall. 

The judgement of the state of convergence from the variation of the l 2 norm 
is more critical in the case of a laminar flow than in the case of a turbulent one. 
That, because in a laminar flow the natural dissipation is much smaller; so when 
the numerical dissipation is reduced at the final stages of a calculation, the overall 
dissipation is small. Thus, instabilities may develop near the surface of a body, 
where due to the fineness of the grid the dissipation is very small. A demonstration 
of the conditions that may develop is provided in figure 11, where the disturbance 
field is shown about the prolate spheroid in the case of a laminar flow. The flow 
conditions for this calculation were a = 10°, Moo = 0.166 and Rei = 7.7x10®. It is 
seen in figure 11 that in a small region near the sting, at the windward side, strong 
disturbances persist, while in the rest calculation domain the disturbances have 
decayed. The persistance of these disturbances affects the quality of the predicted 
flow, locally. If the scaled by the Jacobian I 2 norm is consulted, this local non- 
convergence is not detected (fig. 12a). In contradiction, the unsealed I 2 norm, 
figure 12b, clearly indicates that the flow hits not converged yet. 

A final illustration of the effect of the Jacobian is given in figure 13, where 
the variation of the I 2 norm in the direction normal to the surface of the prolate 
spheroid is shown, at two positions; one on the windward side of the flow (<f> = 30°) 
and one on the leeward side ( <f> = 150°). This figure corresponds to figure 4b. It it 
seen, in figure 13, that in the case of the unsealed by the Jacobian I 2 norms, their 
values are orders of magnitude higher near the surface than in the far field. On the 
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contrary, the scaled Zj norms have higher values in the far field and not near the 
surface. 

Also, noteworthy, in figure 13 is the fact that the I 2 norms have an order-of- 
magnitude-higher values in the region of the vortices, on the leeward side, than in 
the corresponding windward region. This fact is in agreement with the observed 
variation of the velocity disturbance field shown in figure 4b. This justifies the use 
of the velocity disturbance field for judging the convergence state of a numerical 
solution. 


4.3 Partial-Grid Calculations 

The observed nonuniformity of the numerical perturbation field and the con- 
centration of the disturbances in rather small regions, during most of the calculation 
time, suggest that for the estimation of a flow field it is worth applying what may 
be called a partial-grid calculation technique. According to this technique, during 
a certain number of iteration steps the calculation will be restricted to a region 
which will be equal to a part of the overall numerical domain. This region will be 
characterized by the existence of large numerical disturbances, while outside this 
region the level of disturbances will be very small or zero. In this way a numerical 
solution will be accelerated, because, although the toted number of the iterations 
will be equal, or even greater, to the number required for the convergence of the 
solution by application of the standard technique of the full-grid calculation, the 
equivalent number of full-grid iterations will be smaller. 

There are various possibilities for the application of the technique of the partial- 
grid calculation. In case of a nonsteady flow, and assuming that an asymptotic 
steady solution is available, it is necessary to iterate from time to time between the 
partial-grid solution, which will be applied to the varying part of the flow, and the 
full-grid calculation, so that the disturbances will be free to propagate. 

In case of a steady flow calulation a reasonable approach is to reach first the 
full-grid calculation to such a level that the disturbances will be significant only in 
a small domain; then to apply a partial-grid calculation, until “silence” conditions 
are established in the disturbance region. After this, the full-grid calculations may 
be repeated till the flow will be assumed to have converged by means of a mean- 
value criterion. Alternatively, the full-grid calculations may not be repeated, after 
the partial-grid calculations, if the overall convergence level is considered as being 
satisfactory. 


11 


The technique of the partial-grid calculation was validated by being applied 
in the case of the impinging jet. The partial-grid calculations started at N=500 
of the full-grid solution. Then the partial-grid calculation was applied first to a 
cylindrical region that starts from the exit of the nozzle, it extends to the ground 
and encloses the numerically most active part of the jet. The ratio of the partial 
grid to the full grid was 1:3. After the disturbances were significantly reduced (300 
steps), the radius of the cylindrical region was increased, so that the downstream 
region of disturbances was enclosed within the calculation domain. The enlarged 
partial-grid was about half of the total grid. The calculation was continued in the 
enlarged partial grid until the total number of steps became equal to the number 
required for the full-grid solution to converge, i.e., 1700 steps. The equivalent full- 
grid number of steps was only 1000. Thus, with this partial-grid calculation about 
40% of the computation time was saved. 

Some instances of the numerical disturbance field during the partial-grid calcu- 
lation are shown in figure 14; the starting conditions (fig. 14a), the end of the 1/3 
partial-calculation (fig. 14b), the extension of the calculation to the 1/2 of the grid 
(fig. 14c) and the final conditions (fig. 14d). If the full-grid calculation procedure 
was applied, the disturbance field at the equivalent number of iterations would be 
like the one shown in figure 14e. The positive effect of the partial-grid calculation 
is evident if the equivalent figure 14d and 14e, are compared. In addition, the final 
partial-grid results compared well with the final full-grid results (fig. 6e ). 

The calculated flow characteristics of the partial-grid solution were the same 
as those of the full-grid solution. For a qualitative comparison, the development of 
the mean value of the vector A Q is shown in figure 15, for both the full-grid and 
the partial-grid calculation. The mean A Q, in the case of the partial-grid calcu- 
lations is based on the actual domain of calculation. The agreement between the 
final convergence level of the two calculations and the discontinuous increase of the 
disturbance value, observed in the partial-grid calculation curve when the grid was 
extended, are the major points to be noted in figure 15. Additional calculations, 
not presented in this paper, have indicated that the discontinuous increase of the 
numerical disturbance is smaller when the partied grid is so selected that the dis- 
turbances at its boundary eire small. However, if a sequence of partial grids is used, 
as in the example presented above, the overall gain may be positive in case that the 
boundaries of the partial grid are located within a region of disturbances. The rate 
of decay of the disturbances may be high enough to offset the strong perturbation 
which will appear. 

To conclude, the application of the partial-grid calculation technique in the 
case of the jet in crossflow has demonstrated the potential of this technique in 
accelerating the convergence of a numerical solution of a flow. In the simple config- 
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uration of the jet it was possible to reduce the number of iterations by 40%. In more 
complicated cases the gain will be higher. For example, if the jet was issued from 
a wing and not from a cylinder, the grid points which would be included within 
the jet-disturbance region, where the partial-grid calculation was applied, would be 
much less than half of the total number of grid points. It is expected that this 
technique will be beneficial for all flow conditions that include localized complex- 
flow regions. Also, the study of store separation by application of the multiple 
mesh scheme (ref. 10) is another class of problems, which seem promising for the 
application of the partial-grid technique. In this scheme a complex configuration 
is mapped by using a major grid about the main component of the configuration; 
minor overset meshes are used to map each additional component like a store. It is 
expected that the disturbances that are generated during the advance in time of the 
minor meshes (stores) do not cover the whole major grid. So, in such a calculation, 
it will be sufficient each time to estimate only a part of the flow field about the 
major component. 


5. CONCLUSIONS 

In this work, the spatial variation of the numerical disturbances that are gen- 
erated during a numerical solution of a flow, is examined. By considering the 
generation and decay of the velocity disturbances, it is shown that the disturbance 
field is not uniform. More specifically, the decay of the disturbances, or equivalently 
the convergence of the solution, is fast in regions where the structure of the flow is 
simple, but it takes a large number of iterations for the disturbances to decay in 
complex-flow regions, as in separation and vortical areas and in shear flows. 

The observed nonuniformity of the numerical disturbances makes possible the 
reduction of the calculation time by application of the partial-grid calculation tech- 
nique, in which a part of the calculation procedure is applied in selective subregions, 
where the velocity disturbances are large, and not within the whole grid. In the 
particular case of the jet in ground effect that the partial-grid technique was tested, 
a 40% reduction of the number of iterations was obtained. It is expected that this 
technique will prove beneficial in large-scale calculations such as the flow about 
complete aircraft configurations at high angle of attack, the study of vortices gen- 
erated on the tips of a helicopter rotor and the shock boundary-layer interaction 
observed on the upper surface of a wing. 

A secondary issue that has been addressed in the present paper is the effect 
of the Jacobian of the coordinate transformation on the values of the mean- value 
criteria, such as the I 2 norm, in case that the Navier-Stokes equations are written 
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in a generalized coordinate system. It has been shown that in regions in which the 
grid is fine, such as near solid boundaries, the perturbation quantities, on which 
the estimation of the norms is based become infinitesimally small, because in these 
regions the Jacobian has very large values. Thus, in reality the norms reflect only 
changes that happen at the outer boundaries of the computation domain, where 
the value of the Jacobian approaches unity, and not in the whole flow field. This 
situation is corrected if the residual of the equations are unsealed by the Jacobian 
before the norms are estimated. 

In this paper the velocity disturbance field was used for judging the convergence 
state of a numerical solution. Actually this tool proved very useful. Details that 
otherwise would not be detected, like the existence of small instability regions near 
the surface of a body, were easily identified. This type of plot is highly recommended 
in large-scale computing, where a calculation requires many computer hours and it 
is not done all at once. 
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Figure l.-The calculated flow field about the prolate spheroid, (a) surface 
oilflow simulation, (b) crossflow velocity vectors at x/L=0.50. 



Figure 2.-Results of jet in crossflow, (a) trajectories of particles released from 
the jet. (b) velocity vectors in the streamwise symmetry plane. 
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Figure 3.- Variation of A Q during the calculation of the flow about the prolate 
spheroid. 
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Figure 4.-Prolate spheroid calculations. History of development of the numer- 
ical disturbance field, after the solution was disturbed by reducing the numerical 
dissipation terms. (See page 8 for description). 
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Figure 5.-Prolate spheroid calculations. History of development of the numer- 
ical disturbance field, after the solution was disturbed by changing the turbulence 
model. (See page 8 for description). 
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(e) 

Figure 6.- Jet in crossflow calculations. History of development of the numerical 
disturbance field, after the solution was disturbed by changing the cross-section of 
the jet from circular to elliptic. (See page 8 for description). 
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Figure 7.-Prolate spheroid calculations. Vari- 
ation of the transformation Jacobian of the 
coordinates, across the grid. 



Figure 8.-Jet in crossflow calculations. Vari- 
ation of the transformation Jacobian of the 
coordinates, across the grid. 
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UNORM/JACOBI 



Figure 9. -Prolate spheroid calculations. History of the Zj norm, (a) unsealed, 
(b) scaled by the transformation Jacobian. 
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Figure 11. -The numerical disturbance field during the solution of a laminar 
flow about the prolate spheroid. 



NUMBER OF ITERATIONS 

Figure 12. -Laminar-flow about a prolate-spheroid. History of the Zj norm, (a) 
scaled, (b) unsealed by the transformation Jacobian. 






Figure 13.-Prolate-spheroid calculation. Spatial variation of the It norm. 
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(c) 


(d) 



(e) 


Figure 14. -Partial-grid calculation.History of development of the numerical dis- 
turbance field. (See page 12 for description). 


26 





Report Documentation Page 


2. Government Accession No. 


rw\SA 

MtUwviI /Wivurii s.mrl 
‘V * 


t. Report No. 

NASATM 100037 

4. Title and Subtitle 


The Spatially Nonuniform Convergence of the Numerical Solution 
of Flows 


3. Recipient's Catalog No. 


5. Report Date 

November 1987 

6. Performing Organization Code 


7. Author(s) 


Argyris G. Panaras 


9. Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035-5000 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


8. Performing Organization Report No. 

A-88016 


10. Work Unit No. 

505-60 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 


15. Supplementary Notes 

Point of Contact: Argyris G. Panaras, Ames Research Center, M/S 258-1, Moffett Field, CA 94035-5000 
( 415 ) 694-4468 or FTS 4644468 
Permanent Address: Agias Elenis 63, Athens, 15772 Greece 


16. Abstract 

The spatial distribution of the numerical disturbances that are generated during a numerical solution 
of a flow, is examined. It is shown that the distribution of the disturbances is not uniform. In regions 
where the structure of a flow is simple, the magnitude of the generated disturbances is small and their 
decay is fast. However in complex-flow regions, as in separation and vortical areas, large-magnitude 
disturbances appear and their decay may be very slow. The observed non uniformity of the numerical 
disturbances makes possible the reduction of the calculation time by application of, what may be called, 
the partial-grid calculation technique, in which a major part of the calculation procedure is applied in 
selective subregions, where the velocity disturbances are large, and not within the whole grid. This tech- 
nique is expected to prove beneficial in large-scale calculations such as the flow about complete aircraft 
configurations at high angle of attack. Also, it has been shown that if the Navier-Stokes equations are 
written in a generalized coordinate system, then in regions in which the grid is fine, such as near solid 
boundaries, the norms become infinitesimally small, because in these regions the Jacobian has very large 
values. Thus, the norms, unless they are unsealed by the Jacobian, reflect only the changes that happen 
at the outer boundaries of the computation domain, where the value of the Jacobian approaches unity, 
and not in the whole flow field. 


17. Key Words (Suggested by Authorlsll 

Computational fluid dynamics 

Numerical solution of Navier-Stokes equations 

Convergence rate 


18. Distribufron Statement 

Unclassified-Unlimited 


Subject Category-64 


19. Security Classif. (of this report) 

Unclassified 


20. Security Classif. (of this page) 

21. No. of pages 

Unclassified 

29 


22. Price 

A03 



